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Abstract 

We discuss several techniques for the evaluation of the generalised Lyapunov exponents which 
characterise the growth of products of random matrices in the large-deviation regime. A Monte 
Carlo algorithm that performs importance sampling using a simple random resampling step is 
proposed as a general-purpose numerical method which is both efficient and easy to implement. 
Alternative techniques complementing this method are presented. These include the computation 
of the generalised Lyapunov exponents by solving numerically an eigenvalue problem, and some 
asymptotic results corresponding to high-order moments of the matrix products. Taken together, 
the techniques discussed in this paper provide a suite of methods which should prove useful for 
the evaluation of the generalised Lyapunov exponents in a broad range of applications. Their 
usefulness is demonstrated on particular products of random matrices arising in the study of scalar 
mixing by complex fluid flows. 
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I. INTRODUCTION 



Products of random matrices arise in many physical models, of disordered media, of wave 
localisation, and of chaotic dynamics in particular. The main quantity of interest is the 
largest Lyapunov exponent, which gives the rate of exponential growth of the products as 
the number N of factors increases to infinity. The free energy of random Ising models, 
for instance, is given by the largest Lyapunov exponent of a product of matrices, as is the 
localisation length of some random Schrodinger operators. We refer the reader to the book 
by Crisanti, Paladin and Vulpiani [1] for a discussion of these and other applications. 

Often, it is necessary to go beyond the almost-sure, infinite- growth of the matrix prod- 
uct captured by the largest Lyapunov exponent, and examine finite- fiuctuations. These 
are characterised by the distribution of the so-called finite-time (or finite- A?") Lyapunov ex- 
ponents, or equivalently by the generalised Lyapunov exponents i{q), which give the growth 
rate of the qth moment of the norm of the matrix product [e.g. 1-4] . At a mathematical level, 
the generalisation involved entails the passage from the (mutliplicative, non-commutative) 
law of large numbers [5-7] to the corresponding theory of large deviations [8, and references 
therein] . 

One area of applications in which multiplicative large deviations and generalised Lya- 
punov exponents played a central role is the transport, mixing and reaction of constituents 
in complex fluid flows. In the last ten years or so, a number of results have related the 
macroscopic dynamics of scalars and fields in fluid flows to the large-deviation statistics of 
the stretching by these flows [see 9, for an early review]. Speciflcally, the generahsed Lya- 
punov exponents associated with the stretching have been found to control the decay rate 
of purely advected passive scalars [9-14], the spatial distribution of reacting scalars [15, 16], 
the reaction rate of fast reactions [17], the distribution of vorticity in certain turbulent flows 
[18, 19], the clustering of inertial particles [20], the magnetic fleld in kinematic dynamo 
models [21], etc. In most of these applications, the complex fluid flows are modelled by 
random processes which either are white in time (Kraichnan-Kazantsev flows), or consist of 
sequences of independent identically distributed (iid) processes (variously termed renewing, 
renovating, or innovating flows). In the latter case, the stretching is controlled by products 
of iid random matrices of the type considered in this paper. 

In many of these applications, it is necessary to evaluate the generalised Lyapunov expo- 
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nents for specific random matrices. Very few exact results are available, however. As is also 
the case for the usual Lyapunov exponent, given in fact by ^'(0), these are essentially limited 
to matrices satisfying an isotropy property that reduces the problem to scalar multiplication 
[1, 22, 23]. Thus approximations to i{q), either perturbative or numerical, need to be ob- 
tained. Crisanti et al. [1] review several techniques including Cook and Derrida's asymptotic 
results for large sparse matrices [24] , the weak-disorder expansion for near-identity matrices, 
the replica trick (applicable when q is even and positive), and the (heuristic) microcanonical 
estimate. Cycle expansions [25, 26] provide yet another technique. However, these tech- 
niques are limited to special ensembles of matrices: the microcanonical and cycle-expansion 
estimates, for instance, are applicable to ensembles drawn from a small number of matrices. 
There is, therefore, a genuine need for numerical techniques that enable the estimation of 
i{q) for a broad range of matrix ensembles. The main aim of the present paper is to develop 
one such numerical technique and to demonstrate its usefulness by applying it to a few 
examples. 

Several of the papers on fluid mixing cited above contain numerical evaluations of the 
generalised Lyapunov exponents corresponding to simple renewing flows, and in particular 
to the alternating sine map [27] that has become a standard tool of the field. Most of 
these estimates are obtained using a straightforward Monte Carlo sampling of either the 
probability distribution of the finite-time Lyapunov exponents, or of the gth moments of 
the norm of the matrix product. This approach, which we refer to as brute-force Monte 
Carlo in what follows, is highly inefficient unless \q\ is small. This is because it attempts 
to sample events that have an exponentially small probability as N ^ oo. Clearly, what is 
needed is some form of importance sampling, which focuses the computational effort on the 
realisations dominating the estimate of i{q). We propose and test a simple algorithms that 
has this property. This algorithm, which we call Resampled Monte Carlo (RMC), falls in 
the category of sequential importance-samphng [28] or 'go-with-the-winners' strategies [29] 
used extensively in statistical physics and elsewhere; it consists of a simple modification of 
the brute-force computation adding a (random) resampling step which drastically reduces 
the sample variance. As a result, it yields accurate estimates of £{q) with ensembles that 
are orders of magnitude smaller than those required for the brute-force estimation. The 
algorithm is very close to the cloning/pruning algorithm recently developed to estimate 
large-deviation statistics for more general Markov chains [30, 31]. However, our focus on 



3 



products of matrices leads to an algorithm that is particularly simple to implement and to 
analyse. 

Recently, Haynes and Vanneste [14] used an alternative approach to the brute-force Monte 
Carlo sampling for the evaluation of i{q) for the alternating sine flow [see also 16]. This ap- 
proach relates i{q) to the eigenvalue of an (infinite-dimensional) eigenvalue problem that can 
be discretised and solved numerically, at least for 2 x 2 and perhaps 3x3 matrices. We review 
this approach, first to compare its results with those of our RMC algorithm, but mostly be- 
cause the eigenvalue problem can be used to derive interesting properties of i{q)- One such 
property relates the function i{q) associated with an ensemble of matrices A to the corre- 
sponding function associated with the complementary ensemble of matrices A^"*^/] det Al^/"?. 
This relationship is of great practical interest since considering dct instead of 

A can lead to more accurate estimates of i{q) for some value of q. We demonstrate the 
usefulness of this observation in some examples. 

For large \q\, i{q) is controlled by exceedingly rare realisations of the matrix products, 
and hence it is difficult to estimate reliably using Monte Carlo numerical methods, even with 
importance sampling. An alternative, which we pursue in this paper, is to take advantage of 
the largeness of |g| to derive asymptotic estimates. Starting with the eigenvalue problem and 
using a WKB ansatz, we obtain the asymptotics of i{q) for ensembles of bounded matrices 
and for matrices with (not necessarily independent) Gaussian entries. These asymptotic 
estimates, together with the RMC method, the eigenvalue formulation, and the replica 
approach (which we briefly discuss) provide a suite of methods which should prove useful 
for the evaluation of the generalised Lyapunov exponents of products of random matrices 
arising in a broad range of applications. 

The plan of the paper is as follows. In section II, we review the definition of the generalised 
Lyapunov exponents £{q) and their connection with the large-deviation distribution of the 
finite-time Lyapunov exponent. We also derive the eigenvalue problem from which £{q) 
can be inferred, and we use it to relate i{q) obtained for the matrices A to its counterpart 
obtained for the matrices det The RMC algorithm is presented and analysed 

in section 111; there we show that the algorithm leads to an unbiased estimate for the 
qth moment of the matrix product, and we derive an expression for the variance of this 
estimate. Section IV is devoted to alternative methods for the evaluation of i{q), namely 
the numerical solution of the eigenvalue problem, the replica method, and the large-j^j 
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asymptotic results. All the methods dicussed in the paper are tested on three examples of 
random-matrix ensembles in section V. The paper concludes with a discussion in section VI. 
A pseudocode implementing the RMC method, and some technical derivations are relegated 
to three appendices. 



II. GENERALISED LYAPUNOV EXPONENTS 
A. Definitions and basic properties 

We consider successive products of a vector Xq e M.'^ by iid random matrices A„ e 
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n — 1,2, ■ ■ • , N. In other words, we consider the recurrence 

X„ = ^„X„_i, n^l,2,---,N. (1) 

We assume that Xq is determistic and normalised: Xq = xq with ||a;o|| = 1. The randomness 
of the matrices implies the choice of a probability measure on M'^^''. We will not be 
specific as to the properties of this measure; what we have in mind, as illustrated by the 
examples of section V, are random matrices defined by a number of random parameters 
taken from smooth distributions such as the normal or uniform distributions. 

Our focus is on the large- A" behaviour of ||Ar;v||- This can be characterized by considering 
the generalised Lyapunov exponents 

hm i-logE||X^,||^ (2) 

Af-j-oo iV 

where E denotes the expectation over the random matrices. Note that these exponents 
are independent of Xq for almost all Xq and realisations of the matrices An [e.g. 1, 4]. 
Correspondly, the large- A^ asymptotics of the moments of H-^^atH is given by 

E||A:^||^~c,e^^(«) (3) 

for some Cg. Note that in the commutative case d = 1, (3) is exact with Cg = 1. An 
alternative to the definition (2) of £{q) that makes the independence on xq obvious is 

£(g)= hm i-logEp^---A||^ (4) 

N^oo iV 

where the matrix norm is the 2-norm, so that \\An ■ ■ ■ Ai\\ is the largest singular value of 
An---Ai. 



The generalised Lyapunov function i{q), sometimes termed free energy, obviously satisfies 
•^(0) = and can be shown to be convex. It is directly related to the statistics of for 
N ^ 1 [e.g. 1, 4]. These are usually described in terms of the (largest) finite- A?" Lyapunov 
exponents 

i^JV-llog||^^^•••A||. (5) 
The large-deviation theory asserts that the pdf Pat of is approximately 

p^(/i)xe-^^W, (6) 

where x denotes rough asymptotic equivalence, that is, asymptotic equivalence of the log- 
arithms as ^ oo. The function g, variously termed rate function, Cramer function or 
entropy, is convex. It attains a minimum at the Lyapunov exponent h, which satisfies 

lim Hn (7) 

for almost all realisations of the random matrices, and it can be taken such that g(h) = 
g'{h) = 0. Note that g is in fact independent of the norm chosen for A]^ ■ ■ ■ Ai, and that 
the same g would be obtained if (5) was replaced by H]si — log ||-^Ar||- Using the latter 
point, Laplace's method can be applied to write 

E||X„|,./e-e-''-)d.xe— 

and conclude from (2) that ^ and g are Legendre transforms of one another, 

l{q)^su^{qh-g{h)). (8) 

h 

(Rigorous conditions on the probability measure for the A„ that guarantee that (6) and (8) 
hold are given in Ref . [8] .) Since g'(h) — Q, the Legendre relationship l'{q) — h gives 

h^e'{0). (9) 

B. Eigenvalue problem 

The generalised Lyapunov exponents i{q) can be found by solving a family of eigenvalue 
problems parameterised by q. To see this, we consider 

ii„(x) = E/(A„---Aix), (10) 



for a given function / : IR*^ — > R. We now derive a backward equation for Un by noting that 

Un+i{x) = E/(^„+i---Ax) = Ef{A^---AiAx) = Eu^{Ax), 

where the last expectation involves the single matrix A only. Thus, for an arbitrary /, the 
Un satisfy the recurrence relation 

Un+i{x) = Eun{Ax), with uo{x) = fix). (11) 

In the particular case where f{x) — so that Un{xQ) — E H-'^nll^, (11) admits solutions 

of the form 

Un{x) = A"||a;||^^;(e), (12) 

where e = e S'^~^ is a unit vector. The scalar A and function v are determined by 

introducing (12) into (11) to obtain 

(V)(e) = A^(e), (13) 
where we have introduced the linear operator Cq defined by 

{Cqv){e) = E\\Aerv{Ae/\\Ae\\). (14) 
Comparing E H-'^rill^ = Unix^) with (2) gives the following: 

Proposition 1 The generalised Lyapunov exponent £{q) is the logarithm of the largest eigen- 
value Ai of (13): 

£(g)=logAi. (15) 

Here we assume that the point of the spectrum with the largest modulus is an eigenvalue, 
Ai. This can be guaranteed under certain assumptions. (See Ref. [8] where the eigenvalue 
problem (13) is studied in order to estabhsh central-hmit and large-deviation results.) Note 
that since jCg maps positive functions to positive function, Ai > 0. 

The characterisation (15) of the generalised Lyapunov exponents is useful for a number 
of purposes. First, it gives a deterministic method for finding i{q) by solving an eigenvalue 
problem, analytically in simple cases and numerically in less simple cases. Second, the 
eigenvalue formulation can be used to examine the convergence of log E ||^iv||^ as A?" — > oo 
and conclude, for instance, that the convergence is typically exponential, with an error 
proportional to |A2/Ai|^, where A2 is the second largest eigenvalue of (13). Third, the 
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eigenvalue formulation makes it possible to establish some useful properties of £{q) which 
we now discuss. 

In Appendix A, we show that the adjoint of Cg is the operator C_g_d,where Cq is de- 
fined as Cq in (13), but with the matrix A replaced by detyl|-^/^. We then have the 
following useful relationships between generalised Lyapunov exponents of the matrices A, 
A-V|det^|V</ and A-\ 

Proposition 2 Let 

^"W = £"L^'°g^ |'!lft"!4;.^.^!)| r(,)= mnJlogE||A-...^r'r (16) 

Then, 

1. i{q)=i{-q-d), 

2. i{q) = i~{—q — d) if the matrices An satisfy det A„ = 1, 

3. £{q) — i{—q — d) if the matrices A„ are symplectic. 
Note that it follows from the first property that 

£(-c^) =logE|detA|-i (17) 

which extends the well-known observation that [6, 32] 

£{-d)^0 if det^ = l. (18) 

The properties in proposition 2 arc established in Appendix A. They are useful in practice: 
because numerical methods for the estimation of the Lyapunov exponents are more accurate 
when |g| is small, estimates for i{q) with q < —d can be obtained efficiently by evaluating 
£{—q — d). Also, the rephca method (described in section IV B), which provides estimates of 
l{q) for q even and positive, can be used for some negative values of q when the proposition 
2 is exploited. 

As a practical tool for the computation of generalised Lyapunov exponents, the eigenvalue 
problem (13) appears limited to small matrices with c? = 2 or c? = 3, because it requires the 
discretisation of an operator acting on functions of d — 1 variables. (See Refs. [14, 16] and 
section IV A below for some implementations with d — 2.) In the next section we describe 
a Monte Carlo method that does not suffer from this limitation. 



8 



III. RESAMPLED MONTE CARLO 



The simplest Monte Carlo method for the estimation of i{q), which we term brute- force 
Monte Carlo, consists in computing the estimator 

K 

K 



7bf _ }_\^\\A{k) ... Ai}^)^\\q 
— 2-^\\^N ^1 ^0|| , 



fe=l 

where the bracketed superscript indexes K independent realisations of the sequences of 
random matrices A„. Clearly, 

so A^~-'^logZ^^ estimates l{q). This method is hopelessly inefficient, however, unless \q\ is 
small. To see why, note that the variance is of is given by 

varZ^^ = -ivar = -i (E WX^f^ - (E \\X^\\yf) ~ "-^ . 

The convexity of l{q) then implies that exp(A'"£(2g)) ^ exp(2A'"£(g)). So the variance of Z^ 
is exponentially large in N, and a number of realisations K » exp[A'"(^(2g') — 2£(g))] is in 
principle necessary for an accurate estimation of (-{q). 

The inefficiency of the brute-force Monte Carlo estimate stems from the fact that for 
finite g, E||Xjv||'^ is dominated by rare realisations which are undersampled unless K is 
exponentially large. To remedy this, we can resample at each iteration so that the dominant 
contributions to E H-'^ivll^ are represented by more realisations; this is the main idea behind 
sequential importance sampling or 'go- with-the- winners' strategies [28, 29]. We describe 
a particularly simple algorithm for such a resampling strategy which we term Resampled 
Monte Carlo (RMC). 



A. Algorithm 

Like the brute-force Monte Carlo, the algorithm relies on iterations and K realisations, 
calculating Xn^ for n = 1, • • • , and k — 1, - ■ ■ ,K. The difference is that the realisations 
are dependent. Rather than Xn , it is convenient to use the corresponding unit vector 



(fc)i 
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Starting with Eq'^ — xo ior k — 1, ■ ■ ■ , K , the algorithm proceeds iteratively with two steps 



at each iteration n: 



1. Draw K random matrices An , and compute 



^?' = Fifei "^"^ «i" = K"e\r (19) 

W^n K-lW 



2. Resample by letting 

Here the Jn'^^ are independent random variables taking values in {1, • • • , K}, with 



= (20) 



0) ^ 
P{ji''=j)-^^ where /3„ = $:af. (21) 
^" fe=i 



The estimate of E ||Xjv||^ is then given by 

1 



:/3i/92---/5iv. (22) 



Note that the resampling step ensures that, at each iteration n, the weight of each re- 
alisation in the estimate of E is the same. Note also that the resampling is tailored 
to a specific value of q. Unlike in the brute-force Monte Carlo, where the same ensemble 
can be used to estimate i{q) for a range of values of q, the RMC approach requires a new 
sampling for each value of q (although it may be possible to use the same sampling for a 
narrow enough range of q). In several applications, though, i{q) is only required for a single 
value of g [e.g. 10, 13, 14]. 

In Appendix B we give a pseudocode for the RMC algorithm. This illustrates the simphc- 
ity of the algorithm, and should be useful for readers wishing to implement it in a specific 
programming language. 



B. Analysis 

To analyse the algorithm further, we note that the NK random matrices An'' involved in 
the computation form K independent paths consisting of the N matrices that are multiplied 

in succession to obtain each E^^\ These paths are 
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where the random variables In\ n — 1, • • •, A?" are determined by k and by the random 
variables Jn^ according to 



j{k) j{k) _ Mn^) j{k) _ Mn-i) 



The factors o;^^ that are computed along the path that yields are then 



a 



•••A XqW^ 



, n = l,---,7V. (23) 



Note that the distribution of the Ir!'^ is that same as that of the \ since the distribution 
of the latter is independent of k] thus, 

Pn 

For a given realisation of the matrices An^ for n = 1, ■ ■ ■ , N and k = 1, - ■ ■ ,K, the probability 
of a particular path 

aUn) aUn-i) 

is then 

pih,---,jNm = ^^^-^t--^, (24) 

P1P2 ■ ■ ■ Pn 

where 

II A(jn) 40l)^ \\q ^ 

cei^")^ ;"^V and (25) 

Here we abuse notation slightly and use the same symbol a to denote, in (23), a random 
variable that depend on both the An^ and the Jn\ and in (25) one that depend only on the 
An^; the same abuse of notation is made for l3. 

To compute the expected value of functions / of produced by the algorithm, we 
note that the corresponding expectation, E ' is a combination of the expectation E over the 
random matrices An^ and of the expectation over the random variables Jn\ Using (24) to 
compute the latter expectation leads to 

E'fiZ^)^E Y: ^/M, (26) 

. , P1P2 ■ ■ ■ Pn 
n,-,jN=i 

where zn — ^1^2 ■ ■ ■ Pn/K^, and the ctn "^ and ^„ defined as in (25). 
Using (26), it is immediate to establish 
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Proposition 3 is an unbiased estimator for E||Xjv||*.' 

E'Z^ = E||X^||''. (27) 

This follows from the computation 

^E\\AN■■■A^xor=E\\XN\\'l, 



which uses (22), (25) and (26). 

In order to estimate the error of Z^v, we obtain the following expression: 



Proposition 4 The expected value of Zj^ is 

h,-,jN=ij[,-,j'N=i 

This is obtained from 



il,---jAr = l 



AT 



on using the definition (25) of the 

Expression (28) makes clear why the variance of Z^ is much smaller than that of Z^. Only 
terms of the K'^^ terms in (28) lead to contributions proportional to exp(A'"^(2g)) (those 
for which ji — j'l, • ' ' 3n — j'n) with all the others leading to much smaller contributions with, 
in particular, {K{K — l))^ proportional to exp(A^£(g)) (those for which ji ^ j[, ■ ■ ■ Jn ^ j'n)- 
In contrast, in E {Z^Y, all the terms arc proportional to exp(A^£(2g)). 

The improvement can be evaluated explicitly in the scalar case d = 1. Admittedly, this is 
an uninteresting case as far as the numerical evaluation of £(g) is concerned, since (2) holds 
exactly for finite N, but it is instructive nonetheless. For d—1, the asymptotic relation (3) 
holds exactly for all N and with Cq—1. It follows that the terms in (28) can be evaluated 
explicitly : if jk = j'^ for / values of k and jk 7^ j'k for the remaining N — I values, 

E ll^lg") • ■■A^^'^xoUA^ji'''^ ■ ■■a[''Ko\\'' = e'^(2,)+2(iv-o%) 
Since there are {^i)K^{K - 1)^-' such terms, (28) becomes 

N 



E'Zl^^Y. (^) - l)^-'e'^(2'')+2(iv-o^(?) ^ ^^m,) + (i^ _ i)e2n.)) 
1=0 ^ ^ 



N 
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The variance is then 



where we have introduced 

7^ = e^(29)-2^te) _ 1. (29) 

Thus the relative variance of is 

- (30) 
(EZ^)' K ' ^ ' 

and the Monte Carlo estimation of E ||^iv||^ by requires only that K ^ N rather than 
K » exp[N{i{2q) — i{q))] as is the case for the brute- force Monte Carlo. This drastic 
gain in computational efficiency is expected to apply also for matrices with d > 1: the 
non-commutativity is likely to modify (30) only through the introduction of an A*"- and 
i^-independent factor on the right-hand side. 

Although we have found that the RMC algorithm performs very well for a broad range 
of random-matrix products, it is useful to have alternatives methods of evaluating £{q) at 
one's disposal. This provides independent checks for the RMC results or, in the case of 
asymptotic approximation for \q\ ^ 1, makes it possible to estimate £{q) when the RMC 
approach becomes unreliable. Such alternative methods are discussed in the next section. 



IV. OTHER ESTIMATES 



A. Solving the eigenvalue problem (13) 

For d = 2 or 3, it is practical to compute £{q) as the largest eigenvalue of the eigenvalue 
problem (13) for functions v on S'^~^. Here we describe an implementation for d — 2. In 
this case, the unit vector e can be parameterised by an angle 9, and v can be expanded in 
a Fourier series, which we write as 

M-l 
n=0 

and truncate at some M. A straighforward discretisation of the eigenvalue problem (13) is 
then obtained by collocation at points 9m = 277111/ M, m — 0, ■ ■ ■ , M — 1. This leads to the 
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generalised matrix eigenvalue problem 

Pv = XQv, (31) 
where v = {vq, ■ ■ ■ , vm-i)^, and the M x M matrices P and Q have entries given by 

Pmn = ^\\Ae{emW^''^^'-^ and Qmn = e'-'-, (32) 
where e{9m) = {cos 9^, sin 9^)"^ and Q{9m) is defined by 

cose(^„) \ _ Ae{9m) 



(33) 



The expectation in the definition of P can be computed using a Monte Carlo approach, and 
the eigenvalue problem solved using a standard technique. 

B. Replica method for positive even q 

A useful method, known as the replica trick [see 1, and reference therein], makes it 
possible to compute i{q) for q positive and even by finding the largest eigenvalue of a 
qd X qd (deterministic) matrix. To see how this can be achieved, observe that the g-fold 
tensor product Xn with itself satisfies 

X^'^ = ^r^n-'i, (34) 

where A^'^ is the g-fold Kronecker product of A^ with itself. Taking the expectation then 
leads to 

EXr^E^^EX^^. (35) 

Therefore 

EX^'^xe^^^l/, (36) 

where Hq is the largest eigenvalue of the qd x qd matrix EA®^, and y e R^'' is the corre- 
sponding eigenvector. Since for g > even, ||^Ar||^ is obtained from by contraction, 

e{q) ^ for q> even. (37) 

The results extends to the case of odd q when the matrices A have only non-negative entries. 
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C. Large- 1 g| asymptotics 

For large l^j, numerical methods that involve taking expectations by sampling become 
inefficient, and it is useful to develop analytic or semi-analytic methods that take advantage 
of |g| » 1 to provide an asymptotic estimate for i{q). The eigenvalue problem (13) is a good 
starting point. Since the expectation is an integral over the random parameters that define 
the matrix ensemble, we can attempt to approximate this integral for |g| ^ 1 using Laplace's 
method. A dominant-balance argument suggests that the eigenfunction v{e), which depends 
implicitly on q, should have the asymptotic WKB form 

v{e) ~ z{e)e'''"^^\ (38) 

where w and z are independent of q. Substituting this into (13) gives 

Az(e)e«"'(^) = Ee«(i°sll^^ll+"'(^^/ll^^ll))^(Ae/||Ae||). (39) 

When the values of ||A|| are bounded, the expectation on the right-hand side is domi- 
nated by the matrices maximising the argument of the exponential (assuming a non-zero 
probability density for the maximising matrices). Concentrating on the case q > 0, this 
gives 

w{e) = sup (log ||Ae|| + w{Ae/\\Ae\\)) - k, (40) 

A 

for some constant k, where the supremum is over the support of the probability measure of 
the random matrices. Note that w is defined up to the addition of an arbitrary constant. 
Equation (40) can be interpreted as a nonlinear eigenvalue problem, with w as the eigen- 
function and K as the eigenvalue. If this eigenvalue problem has a solution, the largest value 
of K governs the rough asymptotics of Ai and hence the asymptotics of i{q), with the result 

e{q) ~ Kq. (41) 

Note that this behaviour implies that the rate function g{h) of the finite-time Lyapunov 
exponents has a vertical asymptote for h — k. Therefore k is also given by the maximum 
possible (largest) finite-time Lyapunov exponent: 

K= lim sup ^log||>liv-A||. (42) 

N-^°°Ai-Am ^ 
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It would of course be difficult to attempt to determine k by sampling the right-hand side 
of this expression. In general, k < sup^ log ||A||, with the equality holding only in special 
cases; see Appendix C. 

The result (41) can be refined by noting that Laplace's method applied to (39) leads to 
the expectation of a term of the form exp(— g(A ^ A^.A — A*)), where A^ is the maximiser 
in (40) and (•, •) is some scalar product (both A^ and (•, •) depend on e). Carrying out the 
expectation yields a factor where D is the dimension of the support of the measure. 

It follows that 

logg + 0(l). (43) 

This asymptotics implies that g ~ —D\og[n — /i)/2, which describes the manner in which 
g{h) approaches the vertical asymptote aX, k — k. 

When ll^ll is unbounded, the matrices A dominating the expectation in (39) are controlled 
by a balance between the argument of the exponential, which grows with and the 

probability density of A which should decrease with \\A\\ if £(g) is to be finite. This means 
that one needs to apply Laplace's method for movable maxima [e.g. 33] and consider the 
g-dependent maximum of 

log||Ae|| +w{Ae/\\Ae\\) + g"Mog7r(A), (44) 

where 7r(A) is the probabihty density of A and w depends on q. For instance, if t^{A) is 
Gaussian, this maximum corresponds to matrices A with 0{q^^'^) entries, leading without 
further calculations to 

^^e^ioggi/^^ hence £(g)~|logg. (45) 
Correspondingly, gr(/i) x e'* for /i » 1. 

V. EXAMPLES 

A. Two-dimensional sine map 

In studies of transport and mixing by complex fiuid fiows, numerous authors have used 
the random sine map proposed by Picrrchumbcrt [27] as a model of a completely chaotic, 
non-divergent flow. In two dimensions, this map is given by 

a:^n+i = a^n + asin(y„-F0i), = y„ 6sin(x„+i ^2), (46) 
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FIG. 1. Generalised Lyapunov exponents of the product of the random matrices (47) with a = b = 
TT. The result of the RMC method (solid line) are compared with those of the brute-force Monte 
Carlo method (dashed line), and of the numerical solution of the eigenvalue problem (dotted line). 
The inset displays a close up of the region — 2 < g < 0. 

where a and b are fixed parameters, and the random angle 01 and 02 are independent 
and uniformly distributed in [0,27r]. The Jacobian matrix d{xn+i,yn+i)/d{xn,yn), whose 
statistics are independent of (a;„,y„), is given at (0,0) by 



It satisfies detA = 1 and hence, since o? = 2, is symplcctic. 

The generalised Lyapunov exponents corresponding to the ensemble of matrices A gener- 
ated by 01 and 02 characterise the separation of nearby particle in the sine fiow. Remarkably, 
their knowledge makes it possible to predict, in some cases at least, the rate of decrease of 
the variance of a passive scalar released in the fiow [10, 11, 14, 19]. Specifically, this rate is 
given by g{0) and, in view of the Lcgcndrc duality of g{h) and £{q), by —i{q*), where is 
such that ^'(g*) = 0. Because of property 3 of proposition 2, q^ = —1. 

In the literature, has been evaluated using brute force Monte Carlo [12, 13] and 

solving the eigenvalue problem [14, 16]. Here we apply the algorithm of section III to 
demonstrate its efficiency. In Figure 1, we compare £{q) obtained for a = 6 = tt using 
different numerical methods: brute force Monte Carlo, RMC, and numerical solution of the 




(47) 
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K = 


500 1000 2000 


N = 20 


0.12 0.052 0.027 


40 


0.24 0.12 0.059 


80 


0.62 0.23 0.10 



TABLE I. Estimation of the normalised variance var ZAr/(E Z^r)^ for the RMC method appUed to 
the matrices (47) for q = —1 and for different values of the number of realisations K and number 
of iterations N. 

eigenvalue problem using 128 Fourier modes and 128 collocation points. For the latter two 
methods, we have used a relatively small ensemble, with K = 1000, while for the brute force 
computation we have used the much larger K = 10^. The number of matrix multiplication 
N was taken as 100 for the RMC but only N = 50 for the brute- force computation which is 
restricted to moderate values of N. Also shown are the very reliable estimates obtained for 
g = 2, 4 and 6 using the replica method. The figure illustrates how impractical the brute 
force computation is to estimate i{q) for, say, q > 2 and g < 0. The other methods, by 
contrast, provide good estimates for a wide range of q. Based on the comparison with the 
replica estimate, the RMC algorithm, which for the parameters chosen is the faster by a 
factor of about 5, appears to be the more accurate method. The inset in the figure zooms 
on the range q e [—2, 0] to emphasise the substantial differences in the estimates in that 
region leading, in particular, the inaccuracy in the estimates of needed for decay-rate 

predictions in the passive-scalar problem. In this regard, we note that a sequence of 500 
RMC computations gives the average and standard deviation = —0.5916 ± 0.0056. 

We have used the example of the two-dimensional sine flow with g = — 1 to assess the 
dependence of the variance of the RMC estimate on the number of realisations K and 
on the number of iterations N. We have estimated this variance by performing 500 com- 
putations of Z]sr for 9 combinations of the parameters K and N. The results are reported 
in Table I. Unsurprisingly, the sample variance scales roughly like 1/K; more interestingly 
it also scales like N in agreement with the behaviour (30) obtained in the scalar case. The 
behaviour (30) can be tested further: since ^(—2) = 0, 7_i = exp(— 2£(— 1)) — 1 ^ 2.26, 
which compares reasonably well with the various estimates of var Zjv/(E Zjv)^ that can be 
obtained from table I. 



18 



25 
20 
15 
10 

5 



2 4 6 8 10 12 

q. 

FIG. 2. Generalised Lyapunov exponents of the product of the random matrices (47) with a = b = tt 
(circles), and a = tt and b = n/S (squares). The result of the RMC method (solid lines) are 
compared with the large-g asymptotic estimate (dashed lines), and the results of the replica method 
(symbols) . 

Returning to figure 1, we note that the estimates of i{q) appear less accurate for negative 
q unless \q\ < 1; this can easily be remedied, however, by using the third property in 
proposition 2, namely £{q) — i{—q — 2), so that the only negative range that needs to be 
considered is q e [—1,0]. 

The estimation of i{q) is truly challenging for large q. Here, we briefly consider it for the 
matrices (47) in order to assess both the reliability of the RMC method, and the asymptotic 
estimate (43). Figure 2 shows i{q) for the matrices (47) with a = b = tt, and with a = ir 
and b = 7r/8. In both cases, i{q) can be approximated according to (43) with D = 2 (since 
the matrices are defined by 2 random angles (pi and 4>2)- The value of k should be derived 
by solving (40). The case a = 6 is special, however. It can be verified in this case that the 
maximum of ||^e|| is achieved for matrices A and unit vectors e such that 74e/||74e|| = e. As 
a consequence, we have that 

K = logsup ||A|| for a = 6. (48) 

A 

This result, which holds for any matrix ensemble such that >le/||Ae|| = e for A and e 
maximizing ||Ae||, is established in Appendix C. It enables a simple evaluation of k when 
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FIG. 3. Iterates w^^\9) for A; = 1,2,3 in the numerical solution of the problem (40) determining 
K for the matrices (47) with a = tt and b = n/8. 

a — h, giving k — 2.467. The corresponding asymptotic estimate (43) is compared in Figure 
2 with the numerical estimates obtained using the RMC and replica methods. The 0(1) 
term in (43) is determined by matching the asymptotic and numerical results for the largest 
value of q on the figure. The figure demonstrates the validity of the asymptotic estimate; it 
also illustrates the reliability of the RMC method (used here with an ensemble of = 1000 
matrices) which provides accurate estimates of £(g) for q as large as 12, at least for matrices 
considered here. 

The simple result (48) is very special. In general, when a^b, the right-hand side of (48) 
is a strict upper bound for k. There is then no explicit expression for and the problem 
(40) must be solved for both k and w{e). We have implemented a numerical solution of this 
problem for the matrices (47). The implementation relies on an iteration: successive iterates 
w^'^^ A; = 1, 2, • • •, regarded as functions of the angle parameterising e, are represented using 
the truncated Fourier series 

M-1 

n=l 

from which the average (n = 0) term is omitted in order to fix the arbitrary constant in the 
definition of w. The iteration scheme 

w^^+^\em) + /^('+'^ = sup [log \\Ae{ern)\\ + wi^^^e^"®^'-) ^ , (49) 
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where Q{9jn) is defined in (33), determines it;('=+^) on the grid points 9jn — 2m7r/M, with 
^(fe+i) f^xed using the condition of zero average for w'^'^'^^\ The supremum is evaluated 
numerically by finding the maximum over a finite number of matrices A obtained for values 
of 01 and 02 on a grid. An inverse Fourier transform then gives wi!!^^^\ and the iteration can 
continue. Figure 3 shows the first three iterates of this method applied in the case a = tt 
and b — n/8. The functions w'^^\9) are defined for 9 — [0, 27r] and vr-periodic; here we show 
an interval of 9 around the maxima of these functions. The first iterate, corresponding to 
the lowest curve, is simply w^^\9) — sup^log ||A(e(^)||. The next two iterates illustrate 
the rapid convergence of the method; after 4 iterations, convergence is achieved, and the 
estimate k = 1.061 is obtained; this is substantially less than logsup^ \\A\\ = 1.385. The 
validity of our asymptotic formula and evaluation of k are confirmed by Figure 2 which shows 
an excellent match between the asymptotic and numerical estimates of i{q)- A similar match 
was found for other values of a and b. 



B. Three-dimensional sine map 



In order to explore matrices that are not symplectic but have determinant 1, we consider 

the stretching by the volume-preserving map of 

Xn+i ^Xn + a sin(y„ + 0i) , Vn+i ^Vn + b sm{zn + 02) , Zn+1 ^ Zn + c sin(x„+i + 03) , (50) 

where the 0j, j = 1, 2, 3 are independent uniformly distributed in [0, 27r]. This map gener- 
alises to three dimensions the two-dimensional alternating sine map (46) . The corresponding 
Jacobian matrix at the origin is 

^ 1 a cos 01 ^ 

1 6 cos 02 • (51) 

y ccos(asin0i -|- 03) ac cos 0i cos(a sin 0i -|- 03) 1 j 

The results of several numerical computations with these matrices are displayed in Figure 
4. In the main panel, we show three different estimates of i{q), all obtained using the RMC 
method with — 100. The first (solid line) applies the RMC method to the matrices A 
with an ensemble size K — 1000; the second (dotted fine) also uses the RMC method but 
with the much smaller ensemble size K — 100. The results illustrate the difficulties that 
arise when evaluating numerically i{q) for q < 0: ior q < —2 in this case, the numerical 
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FIG. 4. Generalised Lyapunov exponents of the product of the random matrices (51) with a = b = 
c = TT. The RMC estimates of i{q) obtained for N = 100 with an ensemble size K = 1000 (solid 
line) and K = 100 (dotted line) are compared with an estimate of l~ {—q — 3) obtained by applying 
the RMC method to with K = 1000 (dashed line). The results of the replica method, applied 
to A (circles) and (squares) are also indicated. The curves in the inset, which displays a close 
up of the region — 3 < g < 0, have been computed using K = 5000. 

estimates appear to be very unreliable, and the situation does not improve much when the 
number of reahsations is increased from K — 100 to K — 1000. The problem is easily 
remedied, however, using property 2 of proposition 2: by applying the RMC algorithm to 
A"^ rather than to A, we estimate i~{q); this estimate, which proves accurate for q > —2, 
then provides a reliable approximation for i{q) with q ^ —1 since i{q) = i~{—q — 3). The 
curve of i~{—q — 3) is shown by the dashed curve in Figure 4. The best estimate of i{q) 
should be read as the dashed curved for g < —1 and the solid curve for q ^ —2. For 
definiteness, one could choose the point q — —d/2 — —3/2 for the transition between the 
two approximations. 

C. Gaussian matrices 

As a last example, we consider the case of Gaussian matrices. When all the entries are 
independent A'"(0, cr^) variables, the statistics of ||^e|| are independent of e, which leads to 
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FIG. 5. Generalised Lyapunov exponents of the product of 8 x 8 symmetric matrices with inde- 
pendent, zero-mean and variance-one Gaussian upper-diagonal entries. The RMC estimates for 
i{q) obtained with K = 2000 (dashed line) and K = 4000 (solid line) realisations arc compared 
with the RMC estimates for i{-q - 8) obtained with K = 2000 (dash-dotted line) and K = 4000 
(dotted line, almost indistinguishable from the dash-dotted line) . The replica estimate for g = 2 is 
indicated by the circle. 

the explicit expression 

£{q) = log E Pe|r = I log(2a2) + log T - log T (^0 , (52) 

for q > —d, with £{q) = oo for q < —d [1]. No such explicit expressions are available when 
the entries are correlated, however, and l{q) needs to be estimated numerically. Here, we 
examine the case of symmetric matrices with iid A'"(0, cr^) upper-diagonal entries. As in 
the case of independent entries, £{q) = oo for g < — d, and so we can expect difficulties 
in estimating l{q) for values of q shghtly larger than —d, say for q < —d/2. It is indeed 
the Figure 5 demonstrates for d = 8: the figure shows the direct estimates for 

i{q) obtained using the RMC algorithm with 2000 and 4000 realisations. The differences 
between the results for q < —d/2 hints at their inaccuracy, as does an examinination of the 
variance of these estimates. More obviously, the estimates fail to capture the rapid growth 
of i{q) as q ^ —d. Once again, we can invoke proposition 2 to remedy this problem, at least 
partially. Applying the RMC algorithm to the matrices A'^/] det A\^/'^ to estimate i{q), then 
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use the equality £{q) — i{—q — d) gives a much better approximation for £{q) in the range 
—d<q'^ —d/2. Thus the approximation for i(q) obtained in this manner with K — 2000 
and K = 4000 are very close to one another and provide a satisfactory estimate for q close 
to g = —d = —8, though the divergence at —d remains difficult to capture. Note that the 
large-g asymptotics (45) has been verified to apply to the symmetric Gaussian matrices with 
d — 8 considered here; it is easy to check directly from (52) that it is satisfied for matrices 
with iid Gaussian entries. 

VI. DISCUSSION 

Motivated by the key role played by the large-deviation statistics of Lagrangian stretching 
in controlling several aspects of fluid mixing, this paper examines the generalised Lyapunov 
exponents of products of independent random matrices. Such products appear in this context 
when the renewing flows, that is, sequences of simple iid steady flows, are used to model 
complex fluid motion. Products of random matrices appear of course in many other areas 
such as disordered media and wave localisation. 

The main aim of the paper is to present and test a reliable numerical procedure for the 
evaluation of the generalised Lyapunov exponents. The procedure proposed remedies the 
undersampling problem that affects the straightforward, brute-force Monte Carlo estimation 
by introducing a resampling step which ensures that the variance of the estimate scales 
linearly with N, the (large) number of matrix multiplications, rather than exponentially. The 
algorithm chosen, which we term Resampled Monte Carlo, is a particularly simple example 
of sequential importance sampling; its efficiency could be improved, e.g. by resampling 
every few iterations only, or by modifying the resampling method [see 28, for alternative 
approaches] . 

In particular, resampling methods can be devised on the model of the PERM method 
used in the simulation of polymer chains [29, 34, 35, and references therein]. In this method, 
the resampling is carried out only for realisations whose weight (i.e., contribution to the 
estimate of E ||X„||'' in our context) exceeds or falls below two chosen thresholds. If a weight 

exceeds the upper threshold, the realisation is cloned a number of times, with the weights 
of each clone divided accordingly; if a weight falls below the lower threshold, the realisation 
is either pruned with probability 1/2 or has its weight doubled. We have implemented a 
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method of this type, using also a random pruning to keep the number of reahsations con- 
stant. The results are similar to those obtained with the RMC method, but the PERM-like 
method proved somewhat slower in the examples we considered. However, we have made no 
attempt at optimising the choice of the parameters that appear in the method (threshold 
values and number of clones). The PERM method has the advantage of potentially alle- 
viating the problem of sample impoverishment which occurs for large N when most of the 
realisations share the same early history. This problem does not appear to be serious for the 
computations of the generalised Lyapunov exponents of the matrix ensembles we treat in this 
paper, because convergence is achieved at moderately large N. Perhaps a more significant 
advantage of the PERM method in our context is that it can be implemented in a depth-first 
version, where the successive matrix multiplications arc performed for a single realisation at 
a time. The drastically reduced memory requirements of depth-first approaches make them 
suitable for the computations of the generalised Lyapunov exponents of very large matrices. 

We have emphasised that the RMC method, and indeed all methods based on 'go-with- 
the- winners' strategies have resampling strategies that are tailored to a particular value of q. 
When estimates of £{q) are desired over a range of values of q, the computational efficiency 
could be improved by using the same ensemble, and hence the same resampling, for several 
values of q within a narrow interval, rather than a separate ensemble for each value of q. 
We do not pursue these improvements here, preferring to focus on the simple version of 
the algorithm which can be easily analysed and already provides a dramatic improvement 
compared with the brute-force Monte Carlo used by many authors. 

In addition to providing a numerical method for the evaluation of the generalised Lya- 
punov exponents, the paper dicusses some of their properties and, in particular, the relation- 
ship between the exponents associated with an ensemble of matrices A and those associated 
with the corresponding ensemble of matrices det This relationship is useful in 

practice to estimate i{q) for negative q, when a direct application of our algorithm to the 
matrices A can be inaccurate. We also examine the asymptotic form of i{q) for \q\ ^ 1 and 
illustrate, on a specific example, how this form can be obtained by semi-analytical means. 
Asymptotic results of this type usefully complement the direct numerical estimates of i{q) 
which require very large samples as |g| increases. 

We conclude this paper by indicating a few possible extensions of the work reported. 
While the paper focuses on the largest generalised Lyapunov exponents, which encode the 
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large-deviation statistics of the largest finite- A?" Lyapunov exponents, analogous statistics 
for lower Lyapunov exponents (proportional to the logarithm of smaller singular values of 
An • • • ^i) are of interest. It would therefore be useful to develop an efficient numerical 
method to evaluate the corresponding generalised Lyapunov exponents. Work along these 
lines is currently in progress. Another useful extension concerns product of correlated ma- 
trices. The algorithm presented in section III uses the independence of the matrices only to 
take a sequential approach, and hence it can be employed for dependent matrices provided 
that the dependence is on the past only, that is, that remain independent of A^ for k > n. 
More involved dependence would require a rethink of the algorithm. Since the literature on 
fluid mixing literature makes extensive use of white-in-time velocity flelds as an alternative 
to renewing flows, it would also be desirable to develop methods for the efficient evaluation 
of generalised Lyapunov exponents in the context of linear stochastic differential equations 
[see 36, for recent analytical results]. Finally, we note that the methods discussed in this 
paper apply to large matrices (d ^ 1) and so could be employed to study the large-deviation 
statistics of discretised infinite-dimensional systems as arise, for instance, in the problem of 
passive scalar decay [37]. 
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Appendix A: Proof of proposition 2 

We first obtain the adjoint in L2{S'^~^) of Cq. Denoting by de the volume element on 
we consider two arbitrary functions v{e) and ^(e) and compute 

/ w{e)Cgv{e)de = E [ \\Ae\\''w{e)v{Ae/\\Ae\\) de 

= E /" \\A-^e'\\-'i-'^w{A-^e' /\\A-^e\\)v{e')de' /\deiAl 
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where we have changed integration variable from e to e' = 74e/||Ae|| and used that de = 
ll^ell'^deVI det A| = ||A-ie'||-<^de7| det A| [cf. 32]. This gives the adjoint of as 

4 = >^-«-<i' (Al) 

where the operator Cq is defined by 

(Cqv) (e) = E {A-'el\\A-'e\\) . (A2) 

Note that Cq is equivalent to £g, with the matrices A replaced by A~^ /\ det Thus, 
according to proposition 1, defined in (16) is the logarithm of the largest eigenvalue of 
Cq. The first part of proposition 2 follows from the fact that Cq and C]^ — C^q^d have the 
same spectrum. The second part is a particular case of the first for det^ = 1. 
We estabhsh the third part of the proposition by showing that 

l{q) — i^{q) for symplectic matrices . 

Recall first that the matrices A are symplectic iff d is even, and 



A^SA = J, where 



I 
-I 



with I the d/2x d/2 identity matrix. Now, we define the operator J acting on functions on 
S'^~^ according to 

iJv) (e) = v{M) 

and compute 

(JCqv) (e) = E ||AJe||% (AJe/||AJe||) = E ||JA-Te||% {U-^e/\\M-^e\\) 

= E \\A-'^e\\''v (M-^e/IIA-^ell) = E \\A-'^e\\'' (Jv) {A-'^e/\\A-'^e\\) 
= {Cfjv) (e), 

where C~'^ is the analogue of Cq with A''^ replacing A. This computation shows that Cq 
and have the same spectrum, hence £(q) = i~'^{q). The results follows from observing 
that i~'^{q) = C-^{q) since they can be expressed as expectation of the largest singular values 
of A'^ ■ ■ ■ A'^^ and ^4^^ • • • A^^ which coincide. 
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Appendix B: Pseudocode for the RMC method 



We give below a pseudocode for the RMC algorithm. The notation is as in section III A 
except for the omission of the superscripts (k) and subscripts n when these are unnecessary 
for the numerical implementation. The variables 7*^^^^ are introduced to perform the random 
resampling (20)-(21) using the uniformly distributed random variables e. 

fix q 

= o,.-.,o),^ = i,-.-,x 

for n = 1 to N 

for k = 1 to K 

draw random matrix A 



(unit vectors in 
(loop over iterations) 
(loop over realisations) 



(fe) 



end 

for A; = 1 to (resampling) 
draw e uniformly in [0, 
j = min{Z e{l,---,K}\ 7^ - e > 0} 

end 

end 

£ = iV-i log {K-^U^^i/3n) (estimate of £{q)) 

end 

Appendix C: Bounds on k 



Starting from (40), we write k as 

K = sup (log \\Ae\\ + w(y4e/||/le||)) - w{e), 

A 

which holds for any e. Taking e — iu,, where maximises w gives 



(CI) 



K < sup log IIAe^^ll, 

A 
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and in particular 

K < sup log ||^e|| = log sup ||^e|| = log sup \\A\\. 

e,A e,A A 

This is also obvious from the fact that Ai < sup^^ II ^^||^. On the other hand, denoting by 
and the maximisers of ||^e||, and evaluating (CI) at e = e*, we have that 

K > log ||Ae*|| + w(74*e*/||A*e*||) - w(e*) = log sup 11^411 + w(^*e*/||74*e*||) - w(e*). 

A 

A consequence is that w(74*e*/||74,,e*||) < w{e^). In the special case where 

||-4*c* II 

the two inequalities obtained imply that 

K, — log sup ||A||. 

A 
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